seurat_dependency (SAVER)

1. Import Data and Library

library(Seurat)

## Attaching SeuratObject

library(SeuratData)

## Registered S3 method overwritten by 'cli':
##   method     from    
##   print.boxx spatstat

## ── Installed datasets ───────────────────────────────────── SeuratData v0.2.1 ──

## ✓ bmcite 0.3.0

## ────────────────────────────────────── Key ─────────────────────────────────────

## ✓ Dataset loaded successfully
## > Dataset built with a newer version of Seurat than installed
## ❓ Unknown version of Seurat installed

library(cowplot)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(SAVER)

# Just an example way to preprocess data, to demonstrate the PCA part
# # You could use SCTransform into RunPCA also
# bm <- suppressWarnings(LoadData(ds = "bmcite"))
# orig_count <- Matrix::t(bm@assays$RNA@counts)
load("CD4_dat.rds")
load("CD4_cor_mat.rds")
# save(cor_mat_list, file = "CD4_cor_mat.rds")
# saver_mat <- Matrix::t(saver_dat$estimate)
# 
# CD4_ind <- which(as.factor(bm@meta.data$celltype.l2) == "CD4 Naive")
# 
# sub_dat <- saver_mat[CD4_ind, ]
# 
# saver_dat <- SAVER::saver(sub_dat)
# save(saver_dat, store_date, high_gene, CD4_ind, file = "CD4_seurat.Rdata")
# store_date <- date()

dim(sub_dat)

## [1] 4500 1884

# first_filter <- apply(sub_dat, 2, function(x) {sd(x)})
# first_filter_ind <- which(first_filter > 0)
# sub_dat <- sub_dat[, first_filter_ind]
# dim(sub_dat)

Apply SAVER

# saver_dat <- SAVER::saver(t(as.matrix(sub_dat)))
# saver_orig <- SAVER::saver(orig_count)
# save(bm, saver_dat, cell_type, sub_dat, sub_celltype, store_date, file = "SAVER_seurat.Rdata")

2-1. Dependency Measures

library(reshape2) # melt function
library(ggplot2) # ggplot function
library(pcaPP) # Fast Kendall function
library(energy) # Distance Correlation
library(Hmisc) # Hoeffding's D measure

## Loading required package: lattice

## Loading required package: survival

## Loading required package: Formula

## 
## Attaching package: 'Hmisc'

## The following objects are masked from 'package:dplyr':
## 
##     src, summarize

## The following object is masked from 'package:SeuratObject':
## 
##     Key

## The following object is masked from 'package:Seurat':
## 
##     Key

## The following objects are masked from 'package:base':
## 
##     format.pval, units

library(zebu) # Normalized Mutual Information
# library(minerva) # Maximum Information Coefficient
library(XICOR) # Chatterjee's Coefficient
# library(dHSIC) # Hilbert Schmidt Independence Criterion
library(VineCopula) # Blomqvist's Beta

make_cormat <- function(dat_mat){
matrix_dat <- matrix(nrow = ncol(dat_mat), ncol = ncol(dat_mat))
  cor_mat_list <- list()
  
  basic_cor <- c("pearson", "spearman")
  # find each of the correlation matrices with Pearson, Spearman, Kendall Correlation Coefficients
  for (i in 1:2){
    print(i)
    cor_mat <- stats::cor(dat_mat, method = basic_cor[i])
    cor_mat[upper.tri(cor_mat, diag = T)] <- NA
    cor_mat_list[[i]] <- cor_mat
    save(cor_mat_list, file = "../temp_cor_data.Rdata")
  }
  
  # functions that take matrix or data.frame as input
  no_loop_function <- c(pcaPP::cor.fk, Hmisc::hoeffd, 
                        VineCopula::BetaMatrix)
  for (i in 3:5){
    print(i)
    fun <- no_loop_function[[i-2]]
    cor_mat <- fun(dat_mat)
    if (i == 4){ # Hoeffding's D
      cor_mat <- cor_mat$D
    }
    cor_mat[upper.tri(cor_mat, diag = T)] <- NA
    cor_mat_list[[i]] <- cor_mat
    save(cor_mat_list, file = "../temp_cor_data.Rdata")
  }
  
  # functions that take two variables as input to calculate correlations.
  need_loop <- c(zebu::lassie, energy::dcor2d, XICOR::calculateXI)

  for (i in 6:8){
    print(i)
    fun <- need_loop[[i-5]]
    
    cor_mat <- matrix(nrow = ncol(dat_mat),
                      ncol = ncol(dat_mat))
    
    for (j in 2:ncol(dat_mat)){
      for (k in 1:(j-1)){
        if (i == 6){
          cor_mat[j, k] <- fun(cbind(dat_mat[, j], dat_mat[, k]), continuous=c(1,2), breaks = 6, measure = "npmi")$global

        } else {
          cor_mat[j, k] <- fun(as.numeric(dat_mat[, j]),
                               as.numeric(dat_mat[, k]))
        }
      }
    }
    
    cor_mat[upper.tri(cor_mat, diag = T)] <- NA
    cor_mat_list[[i]] <- cor_mat
    save(cor_mat_list, file = "../temp_cor_data.Rdata")
  }
  
  return(cor_mat_list)
}

draw_heatmap <- function(cor_mat){
    len <- 5
    melted_cormat <- melt(cor_mat)
    melted_cormat <- melted_cormat[!is.na(melted_cormat$value),]
    break_vec <- round(as.numeric(quantile(melted_cormat$value,
                                           probs = seq(0, 1, length.out = len),
                                           na.rm = T)),
                       4)
    break_vec[1] <- break_vec[1]-1
    break_vec[len] <- break_vec[len]+1
    melted_cormat$value <- cut(melted_cormat$value, breaks = break_vec)
    heatmap_color <- unique(melted_cormat$value)
  
    heatmap <- ggplot(data = melted_cormat, aes(x = Var2, y = Var1, fill = value))+
      geom_tile() +
      ggplot2::scale_fill_manual(breaks = sort(heatmap_color), 
                                 values = rev(scales::viridis_pal(begin = 0, end = 1)
                                              (length(heatmap_color)))) +
      theme_bw() + # make the background white
      theme(panel.border = element_blank(), panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(), axis.ticks = element_blank(),
            # erase tick marks and labels
            axis.text.x = element_blank(), axis.text.y = element_blank())
    
    return (heatmap)
}


make_cor_heatmap <- function(dat_mat, cor_mat_list){
  fun_lable <- c("Pearson's Correlation", "Spearman's Correlation", "Kendall's Correlation",
                 "Hoeffding's D", "Blomqvist's Beta", "NMI",
                  "Distance Correlation", "XI Correlation")
  
  cor_heatmap_list <- list()
  cor_abs_heatmap_list <- list()
  
  # make correlation matrices
  #cor_mat_list <- make_cormat(dat_mat)
  
  for (i in 1:8){
    print(i)
    if ((i == 4) | (i == 5)) {
      cor_heatmap_list[[i]] <- NULL
      cor_abs_heatmap_list[[i]] <- NULL
      next 
    }
    cor_mat <- abs(cor_mat_list[[i]])
    
    # get heatmaps
    cor_heatmap <- draw_heatmap(cor_mat)
    
    # ggplot labels
    ggplot_labs <- labs(title = paste("Heatmap of", fun_lable[i]),
                      x = "",
                      y = "",
                      fill = "Coefficient") # change the title and legend label
      
    cor_heatmap_list[[i]] <- cor_heatmap + ggplot_labs
    
    if (i %in% c(1,2,3,4,6)){
      cor_abs_mat <- abs(cor_mat_list[[i]])
      cor_abs_heatmap <- draw_heatmap(cor_abs_mat)
      ggplot_abs_labs <- labs(title = paste("Abs Heatmap of", fun_lable[i]),
                              x = "", # change the title and legend label
                              y = "", 
                              fill = "Coefficient") 
      cor_abs_heatmap_list[[i]] <- cor_abs_heatmap + ggplot_abs_labs
    } else {
      ggplot_abs_labs <- labs(title = paste("Abs Heatmap of", fun_lable[i]),
                              subtitle = "Equivalent to Non-Abs Heatmap",
                              x = "", # change the title and legend label
                              y = "", 
                              fill = "Coefficient") 
      cor_abs_heatmap_list[[i]] <- cor_heatmap + ggplot_abs_labs
    }
  }
  
  ans <- list(cor_heatmap_list, cor_abs_heatmap_list)
  
  return (ans)
}

store_date <- date()
# save(cormat_list, heatmap_list, saver_mat, store_date,
#     file ="seurat_corr_saver.RData")

# cormat_list <- make_cormat(saver_mat)

load("SAVER_beta_mat.rds")

heatmap_list <- make_cor_heatmap(sub_dat, cor_mat_list)

## [1] 1
## [1] 2
## [1] 3
## [1] 4
## [1] 5
## [1] 6
## [1] 7
## [1] 8

cor_pearson_mat <- cor_mat_list[[1]]; cor_spearman_mat <- cor_mat_list[[2]];
cor_kendall_mat <- cor_mat_list[[3]]; cor_hoeffd_mat <- cor_mat_list[[4]];
cor_blomqvist_mat <- beta_mat; cor_MI_mat <- cor_mat_list[[6]];
cor_dist_mat <- cor_mat_list[[7]]; cor_XI_mat <- cor_mat_list[[8]];

1. Pearson’s correlation coefficient

• Pearson’s correlation is to measure linear dependency of data, X and Y
• \(-1 \leq \rho_{Pearson}(X, Y) \leq 1\)
• \(\rho_{Pearson}(X, Y) = \frac{\sum(x_i-\bar{x})(y_i -\bar{y})}{\sum(x_i-\bar{x})^2(y_i -\bar{y})^2}\)

set.seed(10)
cor_pearson_mat[1:5,1:5]

##               PTMA       RPS5     RPS24    RPL37A RPS17
## PTMA            NA         NA        NA        NA    NA
## RPS5   -0.09727472         NA        NA        NA    NA
## RPS24   0.07422127 0.07593297        NA        NA    NA
## RPL37A -0.15594202 0.63299248 0.4090960        NA    NA
## RPS17  -0.17080615 0.70654795 0.3664607 0.9208068    NA

quantile(cor_pearson_mat, na.rm = T)

##           0%          25%          50%          75%         100% 
## -1.000000000 -0.006735407  0.020630543  0.120213354  1.000000000

quantile(abs(cor_pearson_mat), na.rm = T)

##           0%          25%          50%          75%         100% 
## 9.034230e-09 1.215906e-02 4.459339e-02 1.392134e-01 1.000000e+00

# plot the smallest correlations
cor_pearson_vec <- sort(abs(cor_pearson_mat), decreasing = T)
plot(cor_pearson_vec)

#plot the high correlations
library(hexbin)
library(RColorBrewer)
high_ind <- quantile(cor_pearson_vec, probs = c(0.9), na.rm = T)
high_ind <- cor_pearson_vec[cor_pearson_vec>high_ind]

par(mfrow = c(2,2))
for(i in 1:4){
 idx <- which(abs(cor_pearson_mat) == sample(high_ind, 1), arr.ind = T)
 idx1 <- idx[1,1]; idx2 <- idx[1,2]
 
 bin <- hexbin(sub_dat[,idx1], sub_dat[,idx2], xbins = 40)
 my_colors=colorRampPalette(rev(brewer.pal(11, 'Spectral')))
 plot(bin, colramp = my_colors,
      xlab = paste0(colnames(sub_dat)[idx1], ", (", idx1, ")"),
      ylab = paste0(colnames(sub_dat)[idx2], ", (", idx2, ")"),
      main = paste0("Correlation of ", round(cor_pearson_mat[idx1, idx2], 3)))
}

#plot the lowest correlations
low_ind <- quantile(cor_pearson_vec, probs = c(0.1), na.rm = T)
low_ind <- cor_pearson_vec[cor_pearson_vec<=low_ind]

par(mfrow = c(2,2))
for(i in 1:4){
 idx <- which(abs(cor_pearson_mat) == sample(low_ind, 1), arr.ind = T)
 idx1 <- idx[1,1]; idx2 <- idx[1,2]
 
 bin <- hexbin(sub_dat[,idx1], sub_dat[,idx2], xbins = 40)
 my_colors=colorRampPalette(rev(brewer.pal(11, 'Spectral')))
 plot(bin, colramp = my_colors, legend = F,
      xlab = paste0(colnames(sub_dat)[idx1], ", (", idx1, ")"),
      ylab = paste0(colnames(sub_dat)[idx2], ", (", idx2, ")"),
      main = paste0("Correlation of ", round(cor_pearson_mat[idx1, idx2], 3)))
}

Heatmap

heatmap_list[[1]][[1]]

2. Spearman’s correlation coefficient

• It captures the monotonic relationship between data, X and Y
• \(-1 \leq \rho_{Spearman}(X,Y) \leq 1\)
• \(\rho_{Spearman} = 1 - \frac{6\sum{d_i^2}}{n(n^2-1)}\) where \(d_i\) is the difference between the ranks of \(x_i\) and \(y_i\)

cor_spearman_mat[1:5,1:5]

##               PTMA       RPS5     RPS24    RPL37A RPS17
## PTMA            NA         NA        NA        NA    NA
## RPS5   -0.15266482         NA        NA        NA    NA
## RPS24   0.05467073 0.03729121        NA        NA    NA
## RPL37A -0.19668940 0.56147422 0.3842762        NA    NA
## RPS17  -0.21570545 0.64125769 0.3371530 0.9104563    NA

quantile(cor_spearman_mat, na.rm = T)

##           0%          25%          50%          75%         100% 
## -1.000000000 -0.004696277  0.085982288  0.239263500  1.000000000

quantile(abs(cor_spearman_mat), na.rm = T)

##         0%        25%        50%        75%       100% 
## 0.00000000 0.03577334 0.12441761 0.25469418 1.00000000

# plot the smallest correlations
cor_spearman_vec <- sort(abs(cor_spearman_mat), decreasing = T)
plot(cor_spearman_vec)

#plot the high correlations
high_ind <- quantile(cor_spearman_vec, probs = c(0.9), na.rm = T)
high_ind <- cor_spearman_vec[cor_spearman_vec>high_ind]

par(mfrow = c(2,2))
for(i in 1:4){
 idx <- which(abs(cor_spearman_mat) == sample(high_ind, 1), arr.ind = T)
 idx1 <- idx[1,1]; idx2 <- idx[1,2]
 
 bin <- hexbin(sub_dat[,idx1], sub_dat[,idx2], xbins = 40)
 my_colors=colorRampPalette(rev(brewer.pal(11, 'Spectral')))
 plot(bin, colramp = my_colors, legend = F,
      xlab = paste0(colnames(sub_dat)[idx1], ", (", idx1, ")"),
      ylab = paste0(colnames(sub_dat)[idx2], ", (", idx2, ")"),
      main = paste0("Correlation of ", round(cor_spearman_mat[idx1, idx2], 3)))
}

#plot the lowest correlations
low_ind <- quantile(cor_spearman_vec, probs = c(0.1), na.rm = T)
low_ind <- cor_spearman_vec[cor_spearman_vec<=low_ind]

par(mfrow = c(2,2))
for(i in 1:4){
 idx <- which(abs(cor_spearman_mat) == sample(low_ind, 1), arr.ind = T)
 idx1 <- idx[1,1]; idx2 <- idx[1,2]
 
 bin <- hexbin(sub_dat[,idx1], sub_dat[,idx2], xbins = 40)
 my_colors=colorRampPalette(rev(brewer.pal(11, 'Spectral')))
 plot(bin, colramp = my_colors, legend = F,
      xlab = paste0(colnames(sub_dat)[idx1], ", (", idx1, ")"),
      ylab = paste0(colnames(sub_dat)[idx2], ", (", idx2, ")"),
      main = paste0("Correlation of ", round(cor_spearman_mat[idx1, idx2], 3)))
}

Heatmap

heatmap_list[[1]][[2]]

3. Kendall’s Tau

• It is an alternative method to Spearman’s correlations, identifying monotonic relationships.
• \(-1 \leq \rho_{Kendall}(X,Y) \leq 1\)
• \(\rho_{Kendall}(X,Y) = \frac{\#\;concordant\;pairs - \#\;discordant \;pairs}{0.5n(n-1)}\)

cor_kendall_mat[1:5,1:5]

##               PTMA       RPS5     RPS24    RPL37A RPS17
## PTMA            NA         NA        NA        NA    NA
## RPS5   -0.10267393         NA        NA        NA    NA
## RPS24   0.03655552 0.02490967        NA        NA    NA
## RPL37A -0.13186234 0.39519822 0.2627274        NA    NA
## RPS17  -0.14498545 0.46029968 0.2286474 0.7401223    NA

quantile(cor_kendall_mat, na.rm = T)

##           0%          25%          50%          75%         100% 
## -1.000000000 -0.003906239  0.064871071  0.178346691  1.000000000

quantile(abs(cor_kendall_mat), na.rm = T)

##         0%        25%        50%        75%       100% 
## 0.00000000 0.02806901 0.09346512 0.19069067 1.00000000

# plot the smallest correlations
cor_kendall_vec <- sort(abs(cor_kendall_mat), decreasing = T)
plot(cor_kendall_vec)

#plot the high correlations
high_ind <- quantile(cor_kendall_vec, probs = c(0.9), na.rm = T)
high_ind <- cor_kendall_vec[cor_kendall_vec>high_ind]

par(mfrow = c(2,2))
for(i in 1:4){
 idx <- which(abs(cor_kendall_mat) == sample(high_ind, 1), arr.ind = T)
 idx1 <- idx[1]; idx2 <- idx[2]
 
 bin <- hexbin(sub_dat[,idx1], sub_dat[,idx2], xbins = 40)
 my_colors=colorRampPalette(rev(brewer.pal(11, 'Spectral')))
 plot(bin, colramp = my_colors, legend = F,
      xlab = paste0(colnames(sub_dat)[idx1], ", (", idx1, ")"),
      ylab = paste0(colnames(sub_dat)[idx2], ", (", idx2, ")"),
      main = paste0("Correlation of ", round(cor_kendall_mat[idx1, idx2], 3)))
}

#plot the lowest correlations
low_ind <- quantile(cor_kendall_vec, probs = c(0.1), na.rm = T)
low_ind <- cor_kendall_vec[cor_kendall_vec<=low_ind]

par(mfrow = c(2,2))
for(i in 1:4){
 idx <- which(abs(cor_kendall_mat) == sample(low_ind, 1), arr.ind = T)
 idx1 <- idx[1]; idx2 <- idx[2]
 
 bin <- hexbin(sub_dat[,idx1], sub_dat[,idx2], xbins = 40)
 my_colors=colorRampPalette(rev(brewer.pal(11, 'Spectral')))
 plot(bin, colramp = my_colors, legend = F,
      xlab = paste0(colnames(sub_dat)[idx1], ", (", idx1, ")"),
      ylab = paste0(colnames(sub_dat)[idx2], ", (", idx2, ")"),
      main = paste0("Correlation of ", round(cor_kendall_mat[idx1, idx2], 3)))
}

Heatmap

heatmap_list[[1]][[3]]

4. Hoeffding’s D

• It tests the independence of data by calculating the distance between the product of the marginal distributions under the null hypothesis and the empirical bi-variate distribution.
• \(-1 \leq D(X,Y) \leq 1\)
• \(D(X,Y) = \frac{(n-2)(n-3)D_1+D_2-2(n-2)D_3}{n(n-1)(n-2)(n-3)(n-4)}\)
  • \(D_1 = \sum_{i=1}^{n} Q_i(Q_i-1)\)
  • \(D_2 = \sum_{i=1}^{n} (R_i-1)(R_i-2)(S_j-1)(S_j-2)\)
  • \(D_3 = \sum_{i=1}^{n} (R_i-2)(S_i-2)Q_i\)

5. Blomqvist’s Beta

• It measures dependency between variables by constructing a two-way contingency table with the medians of each margin as cutting points.
• \(0 \leq \beta \leq 1\)
• \(\beta_n = \frac{n_1-n_2}{n_1+n_2} = \frac{2n_1}{n_1+n_2} - 1\)
• \(\beta = P\{(X-\tilde{x})(Y-\tilde{y})>0\} - P\{(X-\tilde{x})(Y-\tilde{y}) < 0\}\)

cor_blomqvist_mat[1:5,1:5]

##           [,1]      [,2]      [,3]      [,4]      [,5]
## [1,] 1.0000000 0.9996054 0.9957819 1.0008901 1.0009946
## [2,] 0.9996054 1.0000000 0.9973966 1.0005082 1.0004887
## [3,] 0.9957819 0.9973966 1.0000000 1.0006723 1.0007602
## [4,] 1.0008901 1.0005082 1.0006723 1.0000000 0.9979941
## [5,] 1.0009946 1.0004887 1.0007602 0.9979941 1.0000000

quantile(cor_blomqvist_mat, na.rm = T)

##           0%          25%          50%          75%         100% 
## -6353.166667  -158.984277   -31.308420    -6.573932  8889.888889

quantile(abs(cor_blomqvist_mat), na.rm = T)

##           0%          25%          50%          75%         100% 
## 9.665521e-06 7.453301e+00 3.322591e+01 1.617147e+02 8.889889e+03

# plot the smallest correlations
cor_blomqvist_vec <- sort(abs(cor_blomqvist_mat), decreasing = T)
plot(cor_blomqvist_vec)

#plot the high correlations
high_ind <- quantile(cor_blomqvist_vec, probs = c(0.9), na.rm = T)
high_ind <- cor_blomqvist_vec[cor_blomqvist_vec>=high_ind]

par(mfrow = c(2,2))
for(i in 1:4){
 idx <- which(abs(cor_blomqvist_mat) == sample(high_ind, 1), arr.ind = T)
 idx1 <- idx[1,1]; idx2 <- idx[1,2]
 
 bin <- hexbin(sub_dat[,idx1], sub_dat[,idx2], xbins = 40)
 my_colors=colorRampPalette(rev(brewer.pal(11, 'Spectral')))
 plot(bin, colramp = my_colors, legend = F,
      xlab = paste0(colnames(sub_dat)[idx1], ", (", idx1, ")"),
      ylab = paste0(colnames(sub_dat)[idx2], ", (", idx2, ")"),
      main = paste0("Correlation of ", round(cor_blomqvist_mat[idx1, idx2], 3)))
}

#plot the lowest correlations
low_ind <- quantile(cor_blomqvist_vec, probs = c(0.1), na.rm = T)
low_ind <- cor_blomqvist_vec[cor_blomqvist_vec<=low_ind]

par(mfrow = c(2,2))
for(i in 1:4){
 idx <- which(abs(cor_blomqvist_mat) == sample(low_ind, 1), arr.ind = T)
 idx1 <- idx[1,1]; idx2 <- idx[1,2]
 
 bin <- hexbin(sub_dat[,idx1], sub_dat[,idx2], xbins = 40)
 my_colors=colorRampPalette(rev(brewer.pal(11, 'Spectral')))
 plot(bin, colramp = my_colors, legend = F,
      xlab = paste0(colnames(sub_dat)[idx1], ", (", idx1, ")"),
      ylab = paste0(colnames(sub_dat)[idx2], ", (", idx2, ")"),
      main = paste0("Correlation of ", round(cor_blomqvist_mat[idx1, idx2], 3)))
}

Heatmap

heatmap_list[[1]][[5]]

## NULL

6. Normalized Mutual Information

• It measures how much one random variable gives information about the other. For example, High mutual information indicates a large reduction in uncertainty.
• \(0 \leq MI(X,Y) \leq 1\), as it is normalized.
• $MI(X,Y) = f_{X,Y} (x,y) log_2 ; dxdy $
• \(MI(X,Y) = \sum \sum p_{X,Y} (x,y) log \frac{p_{X,Y} (x,y)}{P_X(x)P_Y(y)}\)

cor_MI_mat[1:5,1:5]

##            [,1]       [,2]       [,3]      [,4] [,5]
## [1,]         NA         NA         NA        NA   NA
## [2,] 0.07080326         NA         NA        NA   NA
## [3,] 0.02739069 0.05052604         NA        NA   NA
## [4,] 0.07032199 0.19680495 0.09024883        NA   NA
## [5,] 0.07149862 0.27005785 0.09091991 0.3858145   NA

quantile(cor_MI_mat, na.rm = T)

##           0%          25%          50%          75%         100% 
## 0.0002222295 0.0315687819 0.0877562722 0.2224968631 1.0000000000

# plot the smallest correlations
cor_MI_vec <- sort(abs(cor_MI_mat), decreasing = T)
plot(cor_MI_vec)

#plot the high correlations
high_ind <- quantile(cor_MI_vec, probs = c(0.9), na.rm = T)
high_ind <- cor_MI_vec[cor_MI_vec>high_ind]

par(mfrow = c(2,2))
for(i in 1:4){
 idx <- which(abs(cor_MI_mat) == sample(high_ind, 1), arr.ind = T)
 idx1 <- idx[1,1]; idx2 <- idx[1,2]
 
 bin <- hexbin(sub_dat[,idx1], sub_dat[,idx2], xbins = 40)
 my_colors=colorRampPalette(rev(brewer.pal(11, 'Spectral')))
 plot(bin, colramp = my_colors, legend = F,
      xlab = paste0(colnames(sub_dat)[idx1], ", (", idx1, ")"),
      ylab = paste0(colnames(sub_dat)[idx2], ", (", idx2, ")"),
      main = paste0("Correlation of ", round(cor_MI_mat[idx1, idx2], 3)))
}

#plot the lowest correlations
low_ind <- quantile(cor_MI_vec, probs = c(0.1), na.rm = T)
low_ind <- cor_MI_vec[cor_MI_vec<=low_ind]

par(mfrow = c(2,2))
for(i in 1:4){
 idx <- which(abs(cor_MI_mat) == sample(low_ind, 1), arr.ind = T)
 idx1 <- idx[1,1]; idx2 <- idx[1,2]
 
 bin <- hexbin(sub_dat[,idx1], sub_dat[,idx2], xbins = 40)
 my_colors=colorRampPalette(rev(brewer.pal(11, 'Spectral')))
 plot(bin, colramp = my_colors, legend = F,
      xlab = paste0(colnames(sub_dat)[idx1], ", (", idx1, ")"),
      ylab = paste0(colnames(sub_dat)[idx2], ", (", idx2, ")"),
      main = paste0("Correlation of ", round(cor_MI_mat[idx1, idx2], 3)))
}

Heatmap

heatmap_list[[1]][[6]]

7. Distance Correlation

• it is measure to identify non-linear relationships between two random variables with energy distances.
• distance correlation is calculated by dividing the distance covariance between X and Y by the product of their distance standard deviations.
• \(0 \leq dCor \leq 1\)
• \(dCor(X,Y) = \frac{dCov(Y,Y)}{\sqrt{dVar(X)dVar(Y)}}\)
  • \(dCov(X, Y) = \sqrt{\frac{1}{n^2} \sum_{k=1, l=1}^{n} A_{k,l}B_{k,l}}\)
  • \(dVar(X) = dCov(X,X) and dVar(Y) = dCov(Y, Y)\)

cor_dist_mat[1:5,1:5]

##             [,1]        [,2]       [,3]      [,4] [,5]
## [1,]          NA          NA         NA        NA   NA
## [2,] 0.021415666          NA         NA        NA   NA
## [3,] 0.003142149 0.004832672         NA        NA   NA
## [4,] 0.034461122 0.291511184 0.12930026        NA   NA
## [5,] 0.041251873 0.381817622 0.09936682 0.7929289   NA

quantile(cor_dist_mat, na.rm = T)

##          0%         25%         50%         75%        100% 
## 0.000000000 0.001485644 0.010227512 0.041820003 0.911987832

# plot the smallest correlations
cor_dist_vec <- sort(abs(cor_dist_mat), decreasing = T)
plot(cor_dist_vec)

#plot the high correlations
high_ind <- quantile(cor_dist_vec, probs = c(0.9), na.rm = T)
high_ind <- cor_dist_vec[cor_dist_vec>high_ind]

par(mfrow = c(2,2))
for(i in 1:4){
 idx <- which(abs(cor_dist_mat) == sample(high_ind, 1), arr.ind = T)
 idx1 <- idx[1]; idx2 <- idx[2]
 
 bin <- hexbin(sub_dat[,idx1], sub_dat[,idx2], xbins = 40)
 my_colors=colorRampPalette(rev(brewer.pal(11, 'Spectral')))
 plot(bin, colramp = my_colors, legend = F,
      xlab = paste0(colnames(sub_dat)[idx1], ", (", idx1, ")"),
      ylab = paste0(colnames(sub_dat)[idx2], ", (", idx2, ")"),
      main = paste0("Correlation of ", round(cor_dist_mat[idx1, idx2], 3)))
}

#plot the lowest correlations
low_ind <- quantile(cor_dist_vec, probs = c(0.1), na.rm = T)
low_ind <- cor_dist_vec[cor_dist_vec<=low_ind]

par(mfrow = c(2,2))
for(i in 1:4){
 idx <- which(abs(cor_dist_mat) == sample(low_ind, 1), arr.ind = T)
 idx1 <- idx[1,1]; idx2 <- idx[1,2]
 
 bin <- hexbin(sub_dat[,idx1], sub_dat[,idx2], xbins = 40)
 my_colors=colorRampPalette(rev(brewer.pal(11, 'Spectral')))
 plot(bin, colramp = my_colors, legend = F,
      xlab = paste0(colnames(sub_dat)[idx1], ", (", idx1, ")"),
      ylab = paste0(colnames(sub_dat)[idx2], ", (", idx2, ")"),
      main = paste0("Correlation of ", round(cor_dist_mat[idx1, idx2], 3)))
}

Heatmap

heatmap_list[[1]][[7]]

8. Chatterjee’s XI Correlation

• It measures the degree of dependence between the variables with concept of rank.
• \(0 \leq \xi_n \leq 1\)
• \(\xi_n(X,Y) = 1 - \frac{3\sum_{i=1}^{n-1} |r_{i+1} -r_i|}{n^2-1}\)
• \(\xi_n(X,Y) = 1 - \frac{n\sum_{i=1}^{n-1}|r_{i+1}-r_i}{2\sum_{i=1}^{n}l_i(n-l_i)}\)

cor_XI_mat[1:5,1:5]

##             [,1]         [,2]       [,3]      [,4] [,5]
## [1,]          NA           NA         NA        NA   NA
## [2,]  0.01928156           NA         NA        NA   NA
## [3,] -0.01089311 -0.008368139         NA        NA   NA
## [4,]  0.02582162  0.200651395 0.08444043        NA   NA
## [5,]  0.02491406  0.275113559 0.07579901 0.6121926   NA

quantile(cor_XI_mat, na.rm = T)

##           0%          25%          50%          75%         100% 
## -0.054723119  0.002086095  0.015836627  0.048754506  0.994614676

# plot the smallest correlations
cor_XI_vec <- sort(abs(cor_XI_mat), decreasing = T)
plot(cor_XI_vec)

#plot the high correlations
high_ind <- quantile(cor_XI_vec, probs = c(0.9), na.rm = T)
high_ind <- cor_XI_vec[cor_XI_vec>high_ind]

par(mfrow = c(2,2))
for(i in 1:4){
 idx <- which(abs(cor_XI_mat) == sample(high_ind, 1), arr.ind = T)
 idx1 <- idx[1]; idx2 <- idx[2]
 
 bin <- hexbin(sub_dat[,idx1], sub_dat[,idx2], xbins = 40)
 my_colors=colorRampPalette(rev(brewer.pal(11, 'Spectral')))
 plot(bin, colramp = my_colors, legend = F,
      xlab = paste0(colnames(sub_dat)[idx1], ", (", idx1, ")"),
      ylab = paste0(colnames(sub_dat)[idx2], ", (", idx2, ")"),
      main = paste0("Correlation of ", round(cor_XI_mat[idx1, idx2], 3)))
}

#plot the lowest correlations
low_ind <- quantile(cor_XI_vec, probs = c(0.1), na.rm = T)
low_ind <- cor_XI_vec[cor_XI_vec<=low_ind]

par(mfrow = c(2,2))
for(i in 1:4){
 idx <- which(abs(cor_XI_mat) == sample(low_ind, 1), arr.ind = T)
 idx1 <- idx[1]; idx2 <- idx[2]
 
 bin <- hexbin(sub_dat[,idx1], sub_dat[,idx2], xbins = 40)
 my_colors=colorRampPalette(rev(brewer.pal(11, 'Spectral')))
 plot(bin, colramp = my_colors, legend = F,
      xlab = paste0(colnames(sub_dat)[idx1], ", (", idx1, ")"),
      ylab = paste0(colnames(sub_dat)[idx2], ", (", idx2, ")"),
      main = paste0("Correlation of ", round(cor_XI_mat[idx1, idx2], 3)))
}

Heatmap

heatmap_list[[1]][[8]]

quantile_mat <- c()
for (i in 1:length(cor_mat_list)){
  if (i == 4) {
    quantile_mat <- rbind(quantile_mat, c(NA, NA))
    next
    }
  quantile_mat <- rbind(quantile_mat,
                        quantile(abs(cor_mat_list[[i]]), probs = c(0.15, 0.85), na.rm = T))
}
rownames(quantile_mat) <- c("Pearson", "Spearman", "Kendall", "Hoeffding's D",
                            "Blomqvist's Beta","Dist. Corr", "NMI", "XI Corr")

quantile_mat

##                           15%        85%
## Pearson          0.0063573876 0.20926968
## Spearman         0.0175490367 0.33646081
## Kendall          0.0143049742 0.25347375
## Hoeffding's D              NA         NA
## Blomqvist's Beta 0.9893333333 1.00000000
## Dist. Corr       0.0173682245 0.35759836
## NMI              0.0005631287 0.07361977
## XI Corr          0.0031408643 0.08356398

# save(quantile_mat, store_date, file = "CD4_seurat_corrq.RData")